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Abstract 

Methods  to  incorporate  the  effect  of  ambient  currents  in  the  prediction  of  nearshore  wave 
transformation  are  developed.  This  is  accomplished  through  the  construction  of  a  finite-element 
coastal/harbor  wave  model  based  on  an  extended  mild- slope  wave-current  equation  that  includes  wave 
breaking.  Improved  boundary  conditions  are  used  to  provide  more  accurate  forcing  and  to  minimize 
spurious  wave  reflections  from  the  boundaries.  Multiple  nonlinear  mechanisms,  appearing  both  in  the 
governing  equations  and  in  the  boundary  conditions,  are  handled  successfully  and  efficiently  with  iterative 
techniques.  The  methods  are  tested  against  results  from  other  types  of  models  based  on  parabolic 
approximations  or  Boussinesq  equations  for  three  wave-current  problems  of  common  interest  and  varying 
complexity.  While  indicating  good  agreement  in  general,  the  analysis  also  highlights  the  limitations  of 
parabolic  approximation  models  in  case  of  strong  local  currents  and  velocity  shear.  We  also  consider  the 
harbor  engineering  problem  pertaining  to  waves  approaching  an  inlet  with  a  jettied  entrance,  where  wave- 
current  interaction  can  create  a  complex  wave  pattern  that  adversely  affects  small  craft  navigation  and 
causes  scouring.  The  role  of  ebb  and  flood  currents  on  wave  transformation  and  on  breaking  in  the  vicinity 
of  the  inlet  is  investigated  using  the  model  in  conjunction  with  hydraulic  laboratory  data.  It  is  found  that 
although  the  ebb  currents  cause  larger  waves  outside  the  inlet,  much  of  the  wave  energy  is  soon  dissipated 
due  to  breaking;  during  the  flood  tide,  in  contrast,  more  wave  energy  can  penetrate  into  the  inlet  throat. 
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1.  Introduction 

Reliable  prediction  of  wave  motion  in  coastal  areas  is  crucial  to  coastal  engineering 
applications  associated  with  nearshore  morphologic  change  and  harbor/inlet  maintenance. 
Wave  transformation  in  shallow  coastal  areas  occurs  largely  as  a  consequence  of 
bathymetric  variations  and  domain  geometry,  and  several  modeling  methods  that  can 
reliably  predict  the  resulting  refraction,  diffraction,  reflection,  and  dissipation  are  now 
used  in  practice  (e.g.  Vogel  et  al.,  1988;  Kirby  and  Dalrymple,  1994;  Panchang  and 
Demirbilek,  2001).  In  some  areas,  however,  ambient  tidal  and  other  currents  can  be  strong 
and  their  effect  on  wave  transformation  can  be  substantial.  They  create  a  Doppler  shift  and 
cause  wave  refraction,  reflection,  and  breaking,  which  can  result  in  overall  redistribution 
of  wave  energy. 

In  this  paper,  we  develop  techniques  for  simulating  the  effects  of  ambient  currents  on 
wave  transformation  in  coastal  areas  with  arbitrary  geometry  using  the  finite  element 
method.  Our  focus  is  on  phase-resolving  models  that  are  based  on  a  mass-balance 
formulation,  as  opposed  to  the  energy-balance  (phase  averaged)  models  such  as  SWAN 
(Booij  et  al.,  1999;  Ris  et  al.,  1999)  and  STWAVE  (Resio,  1993).  While  energy-balance 
wave  models  are  suited  to  large  scale  wave  growth  and  wave  transformation  applications, 
phase-resolving  wave  models  are  better  suited  to  domains  with  complex  bathymetric  and 
geometric  features,  where  the  effects  of  wave  diffraction  and  reflection  can  be  important. 
Phase  resolving  models  are  based  on  the  elliptic  mild  or  steep- slope  equation  (e.g. 
Berkhoff,  1972;  Chamberlain  and  Porter,  1995),  which  is  applicable  to  the  full  spectrum  of 
water  waves,  or  on  the  Boussinesq  equations  (e.g.  Madsen  and  Sprensen,  1992;  Nwogu, 
1993;  Wei  et  al.,  1995),  which  are  traditionally  limited  to  shallow  water  waves.  Each  of 
these  approaches  is  accompanied  by  its  own  set  of  modeling  difficulties,  as  described  in  the 
review  by  Panchang  et  al.  (1998).  Here,  we  limit  our  approach  to  the  former  type  of 
models,  which  is  based  on  the  extension  of  Berkhoff’ s  mild-slope  equation  (Kirby,  1984) 
that  includes  the  effect  of  ambient  currents. 

Solving  the  elliptic  equation  presents  considerable  computational  difficulties,  as 
detailed  by  Panchang  et  al.  (1991)  and  Grassa  and  Flores  (2001)  involving  wave-current 
interaction  further  complicates  the  problem.  Existing  models  based  on  the  extended  wave- 
current  equation  have,  therefore,  relied  on  various  approximations.  For  example,  Liu 
(1983)  and  Kirby  and  Dalrymple  (1994)  have  used  the  ‘parabolic  approximation’  in  which 
a  solution  can  easily  be  obtained  by  marching  forward  from  one  computational  row  to  the 
next.  However,  this  method  places  restrictions  on  the  direction  of  wave  propagation  and  is 
hence  applicable  to  cases  with  limited  refraction  and  diffraction  and  with  no  reflection. 
Li  and  Anastasiou  (1992)  proposed  solving  the  extended  elliptic  governing  equation  using 
a  combination  of  the  multigrid  finite  difference  method  and  a  logarithmic  transformation 
of  the  velocity  potential.  The  transformation,  as  noted  by  Radder  (1992),  yields  inaccurate 
predictions  near  complex  bathymetries;  and  the  multigrid  method,  like  the  parabolic 
approximation,  is  best  suited  to  modeling  domains  of  simple  shape  such  as  rectangular 
domains.  Further,  if  the  current  pattern  is  complex,  the  effects  of  these  limitations  are 
likely  to  be  exacerbated. 

Here,  we  present  a  comprehensive  treatment  of  several  issues  associated  with  finite- 
element  modeling  of  wave-current  interaction  using  the  mild-slope  equation. 
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(The  treatment  is  identical  for  the  steep-slope  equation  of  Chamberlain  and  Porter,  1995). 
We  tackle  the  full  elliptic  version  of  the  extended  wave-current  equation  with  improved 
treatment  of  boundary  conditions,  intending  to  better  accommodate  two-dimensional 
wave  scattering  phenomena  in  domains  with  arbitrary  geometry  and  complex  current 
fields.  Wave  breaking  due  to  currents  and  bathymetry  is  also  included,  using  the 
formulation  of  Battjes  and  Janssen  (1978). 

The  work  described  here  is  an  extension  of  recent  developments  in  finite-element 
elliptic  wave  modeling  (Kostense  et  al.,  1988;  Xu  et  al.,  1996;  Panchang  et  al.,  2000;  Zhao 
et  al.,  2001;  Chen  et  al.,  2002,  etc.).  The  layout  of  this  paper  is  as  follows.  After  presenting 
the  governing  equations  in  Section  2,  we  develop  appropriate  boundary  conditions  for  two 
kinds  of  problems  in  Section  3.  Section  4  describes  the  solution  method  and  numerical 
techniques  for  handling  multiple  nonlinear  mechanisms.  In  Section  5,  three  test  cases  for 
model  validation  are  described.  In  Section  6,  numerical  results  are  used  in  conjunction 
with  hydraulic  model  data  to  investigate  the  harbor  engineering  problem  involving  waves 
approaching  a  harbor  through  a  jettied  inlet  in  the  presence  of  ebb  and  food  currents. 
Concluding  remarks  are  given  in  Section  7. 


2.  Governing  equations 

Kirby  (1984)  modified  the  earlier  works  of  Booij  (1981)  and  Liu  (1983)  and  derived  the 
following  extension  of  the  mild-slope  equation  that  includes  the  effect  of  background 
currents: 

D2$  -  D0  ?  9 

~dF  +  (V'u)d 1~  v<cc37^  ~  (k  CCz  ~  a  ^  =  0  (D 

where  the  operator  D/D t  =  ( d/dt  +  U  •  P7),  U (. x ,  y)  =  ( u ,  v)  is  the  current  field,  and  &(x,  y,  t ) 
is  the  complex  velocity  potential  at  mean  water  surface.  The  wave  number  k  and  the 
intrinsic  wave  frequency  a  may  be  determined  through  the  linear  dispersion  relation  and 
the  Doppler  shift  relation 

a2  =  gk  tanh (kh)  (2) 

a  =  0)  —  k-U  (3) 

where  h  is  the  water  depth  relative  to  mean  water  level,  to  is  the  apparent  frequency,  k  is 
the  wave  number  vector  taking  into  account  the  direction  of  wave  propagation.  The  phase 
velocity  C  and  the  group  velocity  Cg  are  defined  as 

C  =  a/k;  Cg  =  da/dk  (4) 

In  order  to  better  simulate  nonlinear  waves  in  shallow  waters,  an  alternative  to  (3)  is  the 
amplitude-dependent  dispersion  relation  developed  by  Kirby  and  Dalrymple  (1986) 

a2  =  gk(  1  +  (kA)2F  i  tanh  (kh))  tanh  (kh  +  kAF2 ) 


(5) 
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where  A  is  the  wave  amplitude  and 

_  cosh(4 kh)  +  8—2  tanh 2{kh)  _  f  kh  \4 

8  cosh A(kh)  ’  ^ 2  \sinh(k/0 ) 

For  practical  applications,  wave  breaking  is  an  important  mechanism,  which  may  be 
induced  by  both  shallow  depths  and  strong  opposing  currents.  A  dissipation  term,  is 
therefore,  added  to  the  model  Eq.  (1)  which  yields  (Booij,  1981;  Liu,  1990) 

D20  -  D0  9  9 

— y-  +  (F-C7) - -  V  (CC„V&)  -  (rcce  -  a  +  iaW)&  =  0  (6) 

Dr  D t  B  B 

where  W ,  a  parameter  associated  with  wave  breaking  and  bottom  friction,  is  a  function  of 
the  (unknown)  wave  height. 

For  monochromatic  waves,  the  wave  potential  can  be  expressed  as 

y,  t )  =  <p(x,  y)exp(-iwf)  (7) 

Substituting  (7)  into  (6)  gives 

V  •  (CCg  F0)  —  V-(U(U- V<p))  +  2  i(oU  •  V(j) 

+  ( k2CCg  +m2  -  a2  +  iwV- U  +  iaW)4>  =  0  (8) 

Kostense  et  al.  (1988)  and  Li  and  Anastasiou  (1992)  have  dropped  the  second  term  in 
(8)  under  the  assumption  that  this  term  is  small  compared  to  the  first  one.  Since  solving  this 
term  takes  very  little  extra  computational  effort  in  a  finite  element  approach,  for  generality 
we  have  retained  the  full  version  of  (8)  as  our  governing  equation. 

The  relation  between  surface  elevation  and  wave  potential  may  be  obtained  from  the 
surface  boundary  condition  as 

and  further  simplification  of  (9)  gives  the  following  linear  relation  which  is  more 
commonly  used  (Liu,  1983;  Kirby,  1984), 

ia 

V=~cp  (10) 

8 

Eq.  (10)  assumes  that  wave  amplitude  gradient  (F|0|/|c[)|)  is  small  in  the  direction  of  the 
current  compared  to  the  phase  gradient  F(arg  0).  The  assumption  is  true  only  if  the  back- 
scattered  and  reflected  waves  are  relatively  weak  compared  to  the  incident  waves. 


3.  Boundary  conditions 

Eq.  (8)  is  elliptic  and  requires  boundary  conditions  along  both  physical  boundaries  (i.e. 
shorelines  and  structures)  and  open  boundaries  (i.e.  artificial  water  boundaries).  A  physical 
boundary,  characterized  by  a  pre- specified  reflection  coefficient,  can  either  partially  or 
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(a) 


(b) 


Fig.  1.  Definition  sketch  of  model  domain,  (a)  Open  sea  application;  (b)  coastal/harbor  application. 


fully  reflect  or  fully  absorb  approaching  waves,  while  an  open  boundary  is  always  intended 
to  be  fully  transmitting  to  both  outgoing  and  incoming  waves.  As  shown  in  Fig.  1,  an 
application  may  be  classified  as  either  an  open- sea  application  that  uses  a  full-circle  open 
boundary,  or  as  a  coastal/harbor  application  for  which  a  semicircle  open  boundary  is  more 
appropriate.  Circular  arcs  are  favored  for  open  boundaries  since  the  directions  of  scattered 
waves  are  more  likely  to  be  perpendicular  to  the  boundary.  Detailed  descriptions  of 
boundary  conditions  for  both  application  types  are  presented  below. 


3.1.  Open  sea  problems 


The  wave  potential  0  at  an  open  boundary  may  be  divided  into  a  pre-determined  part 
0o,  and  an  unknown  part  05  that  represents  the  scattered  wave  field.  For  an  open  sea 
problem,  it  is  assumed  that  the  bathymetric  change  and  ambient  currents  outside  the 
modeling  domain  are  small  and  have  no  significant  effect  on  incident  waves.  As  such,  0O 
can  be  specified  as  the  incident  wave  condition  in  a  plane-wave  form 

0O  =  (—ig/a)Ai  exp [ikR  cos(0  —  0^1  (11) 


where  R  is  the  radius  of  the  open  boundary  circle,  ( R ,  0)  denotes  a  node  location  on  the 
boundary  in  polar  coordinates  (Fig.  la),  At  and  0t  are  the  surface  amplitude  and  angle  of  the 
monochromatic  incident  waves.  05  represents  the  effect  of  non-homogeneities  and  is 
induced  by  currents,  bathymetric  variations,  and  physical  boundary  reflections  within  the 
modeling  domain;  to  describe  this,  a  parabolic- type  absorbing  boundary  condition  may  be 
used  (Xu  et  al.,  1996;  Kirby,  1989), 


dxn 


=  p<t>s+q 


dx, 


(12) 
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where  (xn,  xs )  denotes  a  local  coordinate  system  with  xn  normal  (in  the  radial  direction)  and 
tangential  to  the  boundary,  p  and  q  are  two  known  functions  which,  following  Xu  et  al. 
(1996),  may  be  written  in  the  polar  coordinate  system  as 


(13) 


By  combining  the  conditions  for  0O  and  05,  the  open  boundary  condition  for  0  becomes 


dcj) 

dxa 


=  p</>  +  q 


d2(p 

dxj 


+  8 


where 


(14) 


:  ~P<t>o 

dx„ 


'  dx 2 


(14a) 


Any  physical  boundary  that  lies  within  the  modeling  domain  may  be  fully  or  partially 
reflective;  a  general  boundary  condition  may  be  written  as 


d0 

dxn 


^  ik  cos  y 


1  —  Kv  exp (ikfi) 
1  +  Kv  exp (ikfi) 


1  AV,  , 

Ai  dxn) 


(15) 


where  Kr  is  the  reflection  coefficient,  f3  is  the  phase  shift  between  incident  and  reflected 
waves  at  the  boundary,  y  is  the  incident  wave  angle  relative  to  the  normal  direction  (xn)  of 
the  coastal  boundary,  and  Ax  is  the  amplitude  of  the  incident  waves  approaching  the 
boundary.  The  first  term  of  (15)  was  initially  derived  by  Isaacson  and  Qu  (1990),  taking 
into  account  the  partial  boundary  reflection  for  waves  incident  at  an  arbitrary  angle;  the 
second  term,  derived  by  Chen  et  al.  (2002),  accounts  for  the  wave  height  gradient  at  the 
boundary.  The  second  term  may  be  estimated  from  the  linear  wave  ray  theory  (energy 
conservation)  as 


1  dAl  ^  _  1  dO  _  1  dCg 
Ai  dxn  2  dxs  2Cg  dxn 


(15a) 


where  as  shown  in  Fig.  2,  6  is  the  approaching  wave  direction  accounting  for  both  wave 
incident  angle  y  and  boundary  orientation  angle  0,  such  that  6  =  0  +  y .  The  two  terms  on 
the  right-hand  side  of  (15a)  represent,  respectively,  the  effects  of  wave  convergence  or 
divergence  and  bathymetric  change.  It  is  noted  that  (15a)  does  not  account  for  the  energy 
dissipation,  and  so  it  is  invalid  if  waves  were  to  break  before  reaching  the  boundary 
(shoreline).  In  that  event,  though,  (15a)  becomes  insignificant  and  can  simply  be 
neglected. 

Eq.  (15)  is  nonlinear  in  the  sense  that  the  wave  angle  y  is  a  function  of  the  unknown  0. 
Most  conventional  wave  models  have  in  practice  neglected  the  outgoing  wave  angle  by 
setting  y  =  0.  This  often  results  in  spurious  effects,  while  simulating  wave  transformation 
in  harbors  or  jettied  inlets  (Isaacson  and  Qu,  1990;  Beltrami  et  al.,  2001;  Chen  et  al.,  2002). 
Further,  as  demonstrated  by  Chen  et  al.  (2002),  (15a)  eliminates  spurious  reflections 
arising  from  boundary  curvature  and  the  depth  cut-off  at  shoreline  boundaries  (since 
zero  depth  is  not  allowed).  The  boundary  condition  (15)  can,  therefore,  overcome 
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d0 


Boundary 


Fig.  2.  Outgoing  wave  rays  at  a  boundary. 


the  limitations  of  more  approximate  forms  traditionally  used  (e.g.  see  Panchang  and 
Demirbilek  (2001)  for  a  review). 

3.2.  Coastal  and  harbor  problems 

The  semi-circular  open  boundary  used  for  coastal  and  harbor  problems  (Fig.  lb)  is 
comparatively  more  difficult  to  deal  with,  primarily  because  of  nearshore  bathymetric 
change  and  ambient  currents  in  the  exterior  of  the  semicircular  boundary.  In  this  case,  the 
wave  potential  </>0>  which  may  include  both  incident  waves  and  reflected  waves  in  the 
exterior,  will  be  affected  by  the  exterior  bathymetry  and  currents,  it  therefore,  no  longer 
remains  in  the  form  of  simple  plane  waves. 

The  difficulties  resulting  from  the  exterior  bathymetry  can  be  partly  overcome  by 
introducing  a  one-dimensional  (referred  to  as  ID,  hereafter)  model  for  </>0  to  properly 
account  for  the  effect  of  the  exterior  bathymetry  (Panchang  et  al.,  2000).  The  ID  model 
rests  on  the  assumption  that  exterior  depth  varies  only  in  the  cross- shore  direction 
(Fig.  lb).  Together  with  the  approximation  of  k  =  F(arg  </>0)  which  leads  to  Snell’s  Law 
(ky  =  k o  sin  6i  =  constant),  we  obtain  the  following  ID  form  of  the  mild-slope  equation  and 
the  boundary  conditions 


(16) 


x  =  0  at  offshore 


(17) 


v  at  shoreline 


(18) 
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where  the  second  term  in  (18)  may  be  neglected  if  wave  breaking  occurs,  before  waves 
reach  the  boundary,  otherwise  it  is  given  by 


1  dAl 
Ai  dx 


j_aq 

2Cg  dx 


In  Eqs.  (16)— (1 8),  the  variable  0  is  related  to  0O  by 
</>o  =  0Wexp  (ikyy) 


(18a) 


(19) 


where  (x,  y)  refers  to  a  coordinate  system  shown  in  Fig.  lb,  with  the  incident  wave 
condition  specified  at  the  origin,  the  offshore  end  of  the  ID  transect. 

In  practical  applications,  two  ID  solutions  0i  and  02  may  be  sought  to  account  for  the 
bathymetry  differences  between  two  sides  of  the  semicircular  open  boundary.  To  obtain  0O 
along  the  boundary,  the  two  solutions  are  further  interpolated  and  mapped  onto  the 
boundary  as 


00  — 


1  +y/R 


01  to  + 


1  -y/R 


02 to  exp(ifc y) 


(20) 


where  R  is  the  radius  of  the  semicircular  boundary. 

While  the  above  treatment  has  been  shown  to  effectively  handle  the  problem  of  the 
bathymetric  variations  in  the  exterior  (Panchang  et  al.,  2000;  Zhao  et  al.,  2000;  Panchang 
and  Demirbilek,  2001),  the  problem  of  varying  ambient  currents  in  the  exterior  is  more 
challenging  and  has  no  standard  solution  yet.  If  the  current  field  were  to  vary  only  in  the  v 
(cross-shore)  direction  in  the  exterior,  i.e.  U  =  ( u(x ),  v(v)),  its  effect  may  be  incorporated 
into  the  ID  model  for  the  solution  of  </>0-  However,  such  a  current  field  is  rare  in  practice 
and  we  perforce  assume  that  the  effect  of  ambient  currents  on  incident  waves  is  negligible 
in  the  exterior  region.  Diligence  should  be  exercised  in  selecting  the  model  domain  and  in 
specifying  0O. 

For  the  scattered  waves,  depth  variations  that  occur  along  the  semicircular  boundary 
may  render  the  parabolic  condition  (12)  inappropriate  (Panchang  et  al.,  2000);  a  simple 
expression  for  q  as  in  (13)  cannot  be  obtained.  In  such  a  case,  the  lower-order  absorbing 
boundary  condition  (15)  that  incorporates  the  incident  angle  would  be  more  applicable. 
Because  the  open  boundary  is  fully  transmitting  (^fr  =  0),  05  along  this  boundary  may  be 
written  as 


d(ps  f  dO 

—  =  [ik  cos  y  -  — — 
OX„  \  2ox, 


1  dCg 
2Cg  dxn 


(ps 


(21) 


It  is  noted  that,  in  cases,  where  the  bottom  slope  is  mild  and  scatterers  are  weak  and 
originate  from  near  the  center  of  a  large  semicircle,  (21)  can  be  effectively  simplified  as 


dx„  \  2 R  rs 


(22) 


which  is  the  same  form  used  by  Panchang  et  al.  (2000). 
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The  overall  open  boundary  condition  for  0  is  then 


dcj) 

dxn 


1  \  d</>0 

*"2 


(23) 


where  quantities  related  to  0O  are  obtained  from  (19)  or  (20). 

The  coastal  boundary  condition  for  coastal/harbor  problems  is  the  same  as  that  for  open 
sea  problems,  and  is  given  by  (15). 


4.  Numerical  implementation 

Eq.  (8)  together  with  the  dispersion  relations  and  the  appropriate  boundary  conditions 
as  described  above  forms  a  formidable  boundary  value  problem  from  the  viewpoint  of 
numerical  implementation.  Nonlinearities  are  encountered  in  the  governing  equations 
(wave-current  interaction  and  wave  breaking)  and  in  the  boundary  conditions  (when  the 
direction-related  expressions  are  employed).  Nonlinear  mechanisms  can  only  be  treated  in 
an  iterative  manner,  allowing  the  boundary  value  problem  to  be  linearized  and  solved 
during  each  iteration  by  means  of  the  finite-element  method.  Three  aspects  of  the 
numerical  implementation  are  discussed  below. 

4.1.  Determination  of  nonlinear  parameters 


4.1.1.  Wave  direction  and  Doppler  shift 

The  effect  of  current-induced  Doppler  shift  is  incorporated  into  the  kinematic 
parameters  k  and  a  through  the  combined  dispersion  relations  (2)  and  (3)  or  (5)  and  (3).  To 
determine  the  parameters  k  and  cr,  it  is  necessary  to  first  calculate  the  wave  direction.  Wave 
direction  is  traditionally  obtained  from  the  phase  gradient,  which  is 

F(arg  0)  =  F(Im(ln  0))  =  Im(F0/0)  (24) 


therefore,  the  wave  direction  6  is  determined  by, 
d(arg  <p )  /  <3(arg  <p)  ( 1  d<p 

™  9  =  ~dT~ /  =  In'U  Hy 


Iml'i* 

v0  dx 


(25) 


Taking  6  as  the  direction  of  vector  k ,  the  Doppler  shift  relation  (3)  can  then  be 
re-written  as 


a  =  cj  —  k(u  cos  6  +  v  sin  0)  (26) 

Once  6  is  known  the  intrinsic  frequency  a  and  the  Doppler  shifted  wave  number  k  may 
be  obtained  by  solving  the  combined  nonlinear  Eqs.  (2)  and  (26),  or  (5)  and  (26)  if  one  uses 
the  nonlinear  dispersion  relation. 

It  is  noted  that  wave  direction  can  be  defined  only  in  a  wave  field  that  is  mainly 
progressive.  In  some  harbor  areas,  the  wave  field  0  may  contain  multi-directional  wave 
components  due  to  wave  scattering  and  reflection.  The  wave  direction  obtained  from  (25) 
may  then  be  viewed  as  the  primary  direction  of  wave  propagation.  Since  the  parameters 
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k  and  a  are  calculated  on  the  basis  of  this  wave  direction,  the  effect  of  currents  (Doppler 
shift)  on  scattered  waves  moving  in  other  directions  cannot  be  correctly  accounted  for  by 
using  these  parameters.  Strictly  speaking,  in  a  combined  wave-current  environment,  wave 
kinematics  and  the  flowfield  dynamics  are  coupled  and  cannot  be  solved  separately.  If  the 
wave  field  is  dominated  by  waves  that  have  readily  identifiable  directions,  the  effect  of 
currents  on  other  (weak)  wave  components  would  remain  small  so  as  not  to  significantly 
affect  the  accuracy  of  modeling  results.  This  is  true  because  the  scattered  field  is  generally 
weak  when  there  is  no  significant  wave  reflection  from  coastal  boundaries.  In  case  of 
strong  coastal  reflection,  we  have  to  assume  that  reflected  waves  do  not  ride  through  the 
major  current  field  in  directions  largely  collinear  with  flow  directions.  Otherwise  the 
reflected  wave  field  must  be  treated  separately  from  the  forward-propagating  wave  field. 
This  aspect  remains  a  topic  of  further  research. 

4.1.2.  Outgoing  wave  angle  at  a  boundary 

The  absorbing  boundary  conditions  (15)  and  (21)  require  knowledge  of  y,  the  incident 
wave  angle  relative  to  the  normal  direction  at  the  boundary.  For  scattered  waves  at  the 
offshore  open  boundary,  this  angle  can  be  determined  in  a  manner  similar  to  (25)  with 


(27) 


At  a  coastal  boundary,  however,  as  noted  in  Section  3,  approaching  waves  may  be 
partially  reflected.  In  the  presence  of  reflected  waves,  (25)  cannot  be  used  to  find  y  because 
of  the  multi-directionality  in  the  wave  potential  0.  Considering  the  combination  of 
incident  and  reflected  waves,  Chen  et  al.  (2002)  and  Steward  and  Panchang  (2001)  have 
developed  the  following  equation  for  determining  the  approaching  wave  direction  and 
demonstrated  its  robustness  for  use  in  nonlinear  calculations, 


(l-^r)  d(arg  </>)  I  d(arg  <j>) 

1  +  2 Kr  cos (k(l)  +  K]  dxs  /  dxn 


/ 


(28) 


Note  that  (28)  reverts  to  (27)  when  boundary  has  zero  reflection  (^  =  0),  and  (27) 
reverts  to  (25)  if  x  and  y  were  used  instead  of  xn  and  xs.  (28)  may  present  difficulties  when 
0(arg  (j))/dxn  and  0(arg  (j))/dxs  are  both  nearly  zero.  This  implies  a  situation  with  standing 
waves,  which  eliminates  the  need  for  an  establishing  y  for  use  in  (15). 

4.1.3.  Breaking  parameter  W 

A  number  of  formulations  for  the  wave  breaking  parameter  W  can  be  found  in  the 
literature,  but  few  of  them  have  been  well-tested  in  an  environment  with  strong  wave- 
current  interaction.  In  the  absence  of  wave-current  interaction,  Zhao  et  al.  (2000) 
examined  five  different  parameterizations  for  W  and  concluded  that  the  formulations  of 
Battjes  and  Janssen  (1978)  and  Dally  et  al.  (1985)  are  the  most  robust  for  use  with  the 
mild-slope  wave  equation.  Later,  Zubier  et  al.  (2003)  indicated  that  the  formulation  of 
Battjes  and  Janssen  (1978)  provides  a  better  fit  to  field  data  than  that  of  Dally  et  al.  (1985). 
We  therefore  use,  in  this  work,  the  breaking  criterion  of  Battjes  and  Janssen  (1978),  which 
has  also  been  used  by  Booij  et  al.  (1999)  in  their  energy-balance  wave  model. 
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4.2.  Finite  element  formulation  and  solution  of  linear  systems 

The  flexibility  afforded  by  the  choice  of  finite  elements  not  only  permits  the  best  fit  to 
arbitrarily- shaped  boundaries,  but  also  provides  the  most  efficiency  in  grid  generation. 
Locally  wavelength-dependent  grids  (with  a  typical  resolution  of  15-20  grids  per 
wavelength)  of  varying  size  can  be  used,  hence  keeping  the  number  of  grids  to  an 
optimum.  (The  finite-difference  method,  on  the  other  hand,  requires  an  excessive  number 
of  grids  since  the  size  would  be  limited  based  on  the  shallowest  depths;  this  compounds  the 
computational  difficulties).  We  used  a  graphics  software  package  called  the  Surface  Water 
Modeling  System  (Zundel  et  al.,  1998)  to  construct  a  network  of  triangular  finite  elements. 

With  the  inclusion  of  wave-current  interaction,  the  governing  Eq.  (8)  becomes  more 
complicated,  requiring  in  the  finite-element  formulation  use  of  the  more  general 
variational  principle  (Galerkin  approach),  rather  than  the  Jacobian  functional  method 
used  by  Xu  et  al.  (1996).  For  each  of  the  n  nodes  in  the  discretized  domain,  a  linear 
weighting  function  is  chosen  on  basis  of  local  triangular  meshes.  Integration  of  (8)  in 
accordance  with  the  Galerkin  principle  leads  to  a  large  sparse  linear  system 

[K]  {d>}  =  {f}  (29) 

where  [K]  is  the  nXn  stiffness  matrix,  {f }  the  load  vector,  and  { c|) }  the  solution  vector  of 
unknowns  < pt  (i=  1,2 ,...,n).  Details  of  derivation  may  be  found  in  Chen,  2002). 

As  in  the  case  without  ambient  current,  (29)  forms  an  extremely  large  (n  is  typically  at 
an  order  of  105)  complex-valued  system.  Without  ambient  currents,  Panchang  et  al.  (1991) 
solved  (29)  by  applying  a  pre-conditioned  conjugate  gradient  iterative  scheme  to  the 
normalized  equation.  Several  other  iterative  methods  have  since  been  used  with  varying 
success  (Li,  1994;  Oliveira  and  Anastasiou,  1998).  The  iterative  scheme  of  Li  (1994)  is 
probably  the  most  efficient  algorithm  for  solving  the  symmetric  linear  system  resulting 
from  the  mild- slope  equation  or  its  variants  without  currents.  However,  when  currents  are 
included,  the  stiffness  matrix  [K]  is  no  longer  symmetric  or  Hermitian,  as  a  result  of  the 
additional  term  U-Vf  in  (8).  We  found  that  Li’s  algorithm  then  does  not  lead  to 
convergence.  Even  the  stabilized  biconjugate  method  (BiCGStab,  Van  der  Vorst,  1992), 
originally  developed  for  solving  nonsymmetric  systems,  failed  to  yield  a  converged 
solution  for  some  of  the  test  cases  presented  in  the  next  section.  It  is  concluded  based  on 
our  testing  that,  for  a  nonsymmetric  and  non-Hermitian  complex-valued  linear  system 
associated  with  (8),  the  algorithm  of  Panchang  et  al.  (1991)  is  probably  the  most  robust  and 
guaranteed  to  converge. 


4.3.  Iterative  scheme  for  multiple  nonlinear  mechanisms 

As  noted  earlier,  the  parameters  k,  a,  W  and  y  are  nonlinear  and  have  to  be  solved, 
together  with  0,  by  iteration.  In  an  iterative  scheme,  these  nonlinear  parameters  are  pre¬ 
determined  based  on  the  solution  of  0  from  the  previous  iteration  step,  so  that  a  new  linear 
system  for  0  can  be  formed  based  on  the  finite-element  method  and  solved  as  discussed 
above.  The  first  linear  solution  0  may  be  obtained  by  assuming  zero  currents  (U=  0),  no 
breaking  (W=0),  and  y  =  0.  The  nonlinear  iterations  proceed  until  the  following  stopping 
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criterion  is  met: 

Maximum  H^.,1  -  k/>,l;-i  |/|</>mcl  <  s  (30) 

i=l,n 

where  i  is  the  node  number,  j  the  number  of  iterations,  |0inc|  the  amplitude  of  the  incident 
wave  potential  at  offshore  (x  =  0),  and  e  is  a  pre-specified  tolerance  (typically  e=  10  _3). 
Since  |0inc|  is  a  fixed  reference  value  for  all  nodes,  (30)  actually  defines  the  maximum 
tolerable  absolute  error  (solution  difference  between  two  successive  iterations)  for  the 
entire  domain. 

The  behavior  of  solution  convergence  is  highly  case-dependent,  so  termination  may 
require  some  judgment  in  adjusting  the  tolerance  level  s,  or  even  the  use  of  different 
stopping  criteria.  In  iterative  methods,  a  very  common  approach  is  to  use  a  relative  error 
(either  maximum  or  average)  for  checking  convergence.  However,  we  found  that  using  the 
relative  error  imposed  unnecessarily  high  accuracy  on  smaller  waves  in  the  breaking  zone 
while  trying  to  maintain  the  desirable  accuracy  in  the  more  prevalent  non-breaking  area. 
Due  to  the  dissipation  term,  iterative  solutions  in  the  breaking  zone  are  more  likely  to 
remain  oscillating  with  small  amplitudes,  which  could  prevent  convergence/model 
termination.  Our  analysis  indicated  that  using  absolute  error  in  the  stopping  criterion,  as 
given  in  (30),  helps  faster  convergence  without  sacrificing  the  accuracy  of  solution  in  the 
higher  wave  zones.  It  appears  that  (30)  is  a  robust  stopping  criterion,  capable  of  identifying 
a  reliable  solution  in  general. 

Iterative  methods  have  been  used  for  nonlinear  mechanisms  such  as  wave  breaking  (De 
Girolamo  et  al.,  1988;  Zhao  et  al.,  2000),  wave  direction-dependent  boundary  conditions 
(Isaacson  and  Qu,  1990;  Steward  and  Panchang,  2001;  Beltrami  et  al.,  2001),  and  wave- 
current  interaction  (Kostense  et  al.,  1988).  There  has  been  a  concern  that  multiple 
nonlinear  mechanisms  could  slow  down  or  even  impede  solution  convergence.  Our  study 
shows  that  multiple  nonlinear  parameters  may  be  updated  and  iterated  simultaneously 
without  difficulties  in  achieving  convergence.  This  indicates  that  highly  nonlinear  wave- 
current  interaction  problems  can  be  solved  effectively  in  such  a  manner. 


5.  Model  validation 

In  this  section,  three  cases  of  wave-current  interaction  of  common  interest  are 
simulated  using  the  methods  described  above.  The  absence  of  data  or  analytical  solutions 
for  typical  cases  of  two-dimensional  wave-current  interaction  precludes  model  validation 
in  the  usual  sense;  however,  the  results  are  evaluated  for  overall  consistency  by  comparing 
them  to  the  results  of  other  models  (either  parabolic  or  Boussinesq  time-dependent  finite- 
difference  models). 

5.7.  Vortex  ring 

Vortex  rings,  commonly  seen  in  open  oceans  and  coastal  seas,  have  important  effects 
on  marine  physical  and  biological  processes  (e.g.  Mapp  et  al.,  1985;  Liu  et  al.,  1994).  The 
complex  wave  patterns  resulting  from  wave  and  vortex  interaction  can  be  captured  by 
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airborne  or  satellite  imagery.  Liu  et  al.  (1994)  combined  image  processing  and  simple 
wave  modeling  to  identify  the  existence  and  strengths  of  underlying  vortex  rings.  The 
modeling  procedure  described  here  can  be  used  as  a  more  sophisticated  alternative  in 
future  studies  of  vortical  flows  and  ocean  gyres.  Although  this  case  can  be  easily  tackled 
using  simpler  parabolic  approximation  models,  we  describe  it  here  for  two  reasons:  (a)  to 
validate  the  code  we  developed  before  applying  it  to  other  cases,  and  (b)  to  highlight 
certain  features  of  the  wave  pattern  that  occur  in  more  practical  harbor  engineering 
problems  discussed  later  (Section  6). 

We  consider  the  case  of  wave-vortex  interaction  studied  by  Mapp  et  al.  (1985)  and 
Yoon  and  Liu  (1989),  who  describe  the  current  field  as 

Vr  =  0 


C6  exp 


r<R 


r>R{ 


(31) 


where  Vr  and  Ve  are  the  velocity  components  in  polar  coordinates  with  the  origin  at  the 
center  of  the  ring.  Following  Yoon  and  Liu  (1989),  we  set 

C5  =  0.9  m/s;  C6  =  1.0  m/s;  N  =  2;  Rx  =  343.706  m; 


r2  =  384.881  m;  R3  =  126.830  m. 


(32) 


Such  a  vortex  ring  creates  a  shearing  current  (maximum  Ve  =  1  m/s)  all  around  the 
circle.  Since  the  water  depth  is  constant  (h=  10  m),  we  choose  a  circular  computational 
domain  (of  radius  R  =1800  m)  using  the  boundary  conditions  described  in  Section  3.1. 
This  case  provides  a  perfect  test  for  pure  wave-current  interaction,  in  which  plane  waves 
(with  an  incident  amplitude  of  At=  1  m  and  period  T  =  19.43  s)  meet  currents  at  all  angles. 
As  shown  in  Fig.  3a,  the  vortex  ring  is  centered  at  (x0,  yo)  =  (0,  1000  m).  The  domain  is 
evenly  discretized  into  triangular  meshes  with  about  20  nodes  per  wavelength,  generating 
a  total  of  102,171  nodes.  The  simulation  starts  with  a  plane  wave  solution  and  convergence 


(a)  *  (b) 


X  (km)  X  (km) 


Fig.  3.  Model  results  for  wave-vortex  ring  interaction,  (a)  Amplitude;  (b)  surface  elevation. 
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is  reached  after  only  four  nonlinear  iterations  for  wave-current  interaction.  Wave  breaking 
is  not  activated  in  this  case.  The  total  run  time  for  this  simulation  is  about  15  min  on  a 
2.4  GHz  desktop  PC. 

The  simulated  wave  amplitude  and  a  snapshot  of  the  surface  elevation  are  presented  in 
Fig.  3.  As  can  be  seen,  wave  energy  has  been  considerably  redistributed  due  to  the 
disturbance  created  by  the  vortex  ring.  The  effect  of  vortical  flow  on  wave  propagation  is 
complicated  because  progressive  waves  encounter  shear  currents  at  all  angles.  The  vortex 
ring  not  only  affects  the  area,  where  opposing  (left-hand  side)  and  following  (right-hand 
side)  currents  exist,  but  also  a  large  part  of  the  down- wave  region.  The  combined  effect  of 
wave  refraction,  diffraction,  and  reflection  give  rise  to  modulations  of  wave  amplitude  that 
are  neither  symmetric  nor  regularly  oriented. 

Yoon  and  Liu  (1989)  modeled  this  problem  with  a  finite-difference,  phase-resolving, 
parabolic  approximation.  The  overall  wave  pattern  obtained  by  them  (Fig.  13a  in  Yoon 
and  Liu,  1989)  is  similar  to  ours  (Fig.  3a  and  b).  Fig.  4  provides  a  detailed  quantitative 
comparison  of  the  two  model  results  for  six  lateral  transects.  Although  the  two  models 
agree  well  for  this  test  case,  it  is  noted  that  the  method  presented  here  has  greater 
versatility  because  it  does  not  rely  on  the  parabolic  approximation.  In  fact,  some 
limitations  of  the  parabolic  approximation  may  be  captured  in  Figs.  3  and  4.  For  instance, 
due  to  the  highly  oblique  angle  with  which  waves  approach  the  shearing  current  in  the 
vicinity  of  point  A  in  Fig.  3a,  wave  caustics  are  expected  to  occur,  causing  a  small  band  of 
waves  to  be  reflected  toward  the  right  (McKee,  1974;  Peregrine  and  Smith,  1975).  The 
reflection  is  manifested  in  the  form  of  amplitude  modulations  in  our  results  (e.g.  down- 
wave  and  to  the  right  of  point  A  in  Fig.  3,  and  to  the  right  of  B  and  C  in  Fig.  4).  The 
parabolic  approximation  model  is  not  expected  to  simulate  wave  reflection  or  large-angle 
diffraction,  hence  the  model  results  of  Yoon  and  Liu  (1989)  show  a  less  disturbed  solution 
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Fig.  4.  Simulated  relative  wave  amplitudes  for  wave-vortex  ring  interaction  (thick  line:  present  model;  thin  line: 
Yoon  and  Liu  (1989)). 
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to  the  right  of  B  and  C  in  Fig.  4.  This  type  of  reflection  of  waves  that  are  obliquely  incident 
on  a  shearing  current  is  relevant  to  our  simulations  pertaining  to  the  interaction  of  waves 
with  tidal  currents  near  an  inlet  (see  Section  6.4). 


5.2.  Rip  current  on  sloping  beach 


Wave-induced  rip  currents  play  a  significant  role  in  nearshore  morphodynamics.  Waves 
propagating  over  a  rip  current  system  on  a  sloping  beach  experience  both  bathymetric 
shoaling  and  current-induced  refraction-diffraction.  For  an  idealized  case,  the  interaction 
of  a  rip  current  and  shoaling  waves  was  originally  studied  by  Arthur  (1950)  using  a  linear 
geometric  ray  model.  The  problem  considered  has  a  sloping  beach  with  a  uniform  slope  of 
1/50.  The  current  field  based  on  a  local  coordinate  system  (shown  in  Fig.  5)  is  given  by 


„=0.1442vF(4)F(^); 


V  =  -0.5192 


\l-( y  )21 

f(  y  ) 

\76.2/  _ 

V76.2/  t 

■jc/7.62 


F(r)dr 


(33) 


with 


f(i,=iexp(-£)  o4> 

In  accordance  with  the  boundary  formulation  described  in  Section  3.2,  we  use  a 
semicircular  domain  (of  radius  R  =  240  m  and  containing  over  85,000  nodes),  for  which 
currents  far  away  from  the  center  of  the  domain  have  no  cross-shore  component  and  hence  no 
effect  on  normally  incident  waves.  A  normally  incident  wave  of  unit  amplitude  is  specified  at  a 
distance  of  240  m  offshore.  (For  greater  clarity,  only  a  portion  of  the  complete  semicircular 
domain  is  shown  in  Fig.  5  and  others  like  it).  The  incoming  wave  potential  0O  along  the  open 
boundary  is  obtained  by  solving  the  ID  model  as  described  by  Eqs.  (16)— (18)  and  is 


(a)  {  Ai=  l 


Fig.  5.  Modeled  wave  field  in  the  presence  of  a  rip  current,  (a)  Amplitude;  (b)  phase. 
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Fig.  6.  Simulated  relative  wave  amplitudes  using  two  boundary  conditions. 


substituted  into  (23)  to  obtain  the  open  boundary  condition  for  the  2D  model.  Wave  breaking 
is  not  activated  for  consistency  with  other  model  results  (e.g.  Yoon  and  Liu,  1989). 

To  avoid  the  singularity  of  zero  depth  at  the  coastline,  we  place  the  shoreline  boundary 
at  y  =  1  m,  where  the  cut-off  depth  is  0.02  m;  it  is  specified  as  non-reflective  boundary 
(Kr  =  0,  /?  =  0)  using  generalized  boundary  condition  (15).  Because  of  the  large  shoaling 
effect  near  the  shoreline,  the  second  term  in  (15)  plays  an  important  role.  This  can  be  seen 
in  Fig.  6,  where  results  for  simulations  with  two  different  incident  wave  periods  are  shown. 
The  improved  boundary  condition  (15)  is  effective  in  eliminating  spurious  reflections 
resulting  from  the  cut-off  depth  at  the  shoreline.  These  spurious  reflections  are  seen  quite 
significant  for  longer  waves,  even  with  a  very  mild  bottom  slope  (1/50)  in  this  case. 

The  model  simulation  for  T=  8  s  in  this  case  took  five  nonlinear  iterations  and  about 
10  min.  Results  shown  in  Fig.  5  indicate  that  the  opposing  rip  current  has  modified  the 
wavelengths  and  caused  a  notable  bending  of  phase  contours.  As  a  consequence,  wave 
energy  is  focused  and  trapped  to  form  a  high- wave  zone  along  the  centerline  across  shore; 
this  is  accompanied  by  two  symmetric  shadow  regions  on  the  sides.  The  combined  effects 
of  the  rip  current  and  bathymetric  shoaling  cause  amplitude  amplification  to  a  factor  over 
six  in  the  central  area. 

Our  results  are  compared  in  Fig.  7  to  the  parabolic  approximation  results  of  Yoon  and 
Liu  (1989)  along  four  longshore  and  two  cross-shore  transects,  and  again  fairly  good 
agreement  is  seen  in  general.  However,  there  is  an  obvious  mismatch  along  v  =  0  for 
y  <60  m.  Our  results  show  that  the  amplitude  has  a  peak  value  of  5.5  at  approximately 
45  m  offshore,  and  after  a  small  drop,  it  increases  to  about  10  at  the  boundary.  In  other 
words,  the  non-breaking  amplitude  has  nearly  doubled  in  the  shallowest  25  m  of  the 
domain,  where  opposing  current  is  becoming  negligible.  The  results  of  Yoon  and  Liu 
(1989),  on  the  other  hand,  show  a  significant  drop  from  the  peak  value  of  5  down  to  3.2. 
While  using  a  different  parabolic  approximation,  Liu  (1983)  speculated  that  this  amplitude 
drop  was  related  to  the  lateral  component  of  the  current.  We  investigated  this  with 
the  elliptic  finite-element  model  and  found  that  suppressing  the  lateral  velocity  u  made 
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Fig.  7.  Simulated  relative  wave  amplitudes,  rip  current  case  (thick  line:  present  model;  thin  line:  Yoon  and  Liu, 
1989). 


little  difference.  On  the  other  hand,  it  may  be  noted  that  shoaling  alone  causes  a  similar 
doubling  at  other  cross-shore  transects  (e.g.  x  =  143.34  m).  We  submit,  therefore,  that  the 
increase  in  the  amplification  obtained  by  the  elliptic  model  is  a  reasonable  feature  that  the 
parabolic  approximation  cannot  capture.  Parabolic  approximations  assume  that  the  second 
derivative  of  wave  amplitude  in  the  dominant  wave  propagation  direction  is  of  higher- 
order  and  hence  dropped;  in  this  case,  the  assumption  may  not  be  appropriate.  The 
mismatch  between  the  parabolic  approximations  and  elliptic  equation  solutions  (Fig.  7) 
are  significant  only  near  the  central  nearshore  area  or  around  the  shadow  areas  (x=  ± 
15  m),  where  wave  amplitudes  apparently  experience  considerable  change  in  the  v 
direction.  At  distant  cross-shore  sections,  the  two  models  agree  because  changes  in  the 
v-direction  are  less  enhanced. 


5.3.  Rip  current  on  beach  with  sandbar 

We  consider  a  second,  more  realistic,  rip  current  case  that  has  been  studied  by  Chen 
et  al.  (1999).  Fig.  8  shows  the  bathymetry  of  a  plane  beach  with  a  discontinuous  offshore 
sandbar  and  the  underlying  rip  current  field.  The  current  field,  provided  to  us  by  Professor 
Q.  Chen  of  the  University  of  South  Alabama,  is  more  complicated  than  that  in  Section  5.2. 
Strong  vorticity  characterizes  much  of  the  current  field,  for  instance,  on  both  sides  of 
the  central  jet  of  the  rip  current  and  also  behind  the  sandbar  in  the  surf  zone.  The  vortices 
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Fig.  8.  Plane  beach  with  gapped  sandbar  (following  Chen  et  al.,  1999)  bathymetry  and  the  underlying  rip  current 
(double-dashed  lines  indicate  sandbar  locations). 

have  spatial  scales  comparable  to  the  wavelengths,  indicating  large  velocity  gradients.  In 
this  case,  the  semicircular  computational  domain  (R=  10  m)  contains  75,871  nodes.  A 
normally  incident  wave  train  with  period  T=  Is  and  amplitude  At  =  0.024  m  is  prescribed 
at  10  m  offshore,  and  the  coastal  boundary,  specified  as  fully  absorbing,  is  placed  at  a  cut¬ 
off  depth  of  0.03  m.  Nonlinear  iterations  are  performed  simultaneously  for  all  nonlinear 
mechanisms,  including  wave  breaking.  In  general,  the  number  of  nonlinear  iterations 
required  is  greater  when  the  breaking  mechanism  is  activated.  This  case  required  13 
iterations  to  reach  convergence,  but  the  total  run  time  is  only  about  8.4  min. 

The  results  of  the  simulated  wave  field  are  shown  in  Fig.  9.  The  snapshot  of  the  surface 
wave  elevation  (Fig.  9a)  reveals  a  complex  wave  pattern  due  to  the  combined  effects  of 
varying  bathymetry  and  underlying  rip  current.  This  wave  pattern  is  visually  similar  to  those 
obtained  by  Chen  et  al.  (1999,  Plate  6)  using  a  modified  Boussinesq  model  (Wei  et  al.,  1995; 
Chen  et  al.,  2000;  Kennedy  et  al.,  1998).  The  resulting  amplitude  and  phase  diagrams 
(Fig.  9b  and  c)  are  understandably  complex.  It  is  clear,  though,  that  along  the  centerline 
large- waves  build  up  and  then  break  right  in  front  of  the  gap,  at  a  depth  of  approximately 
0.12  m.  Away  from  the  central  channel,  the  breaking  line  is  largely  along  the  toe  of  the 
underwater  sandbar.  Behind  the  bar,  waves  appear  to  be  disturbed  by  the  vortical  flow, 
which  results  in  fluctuating  phase  lines  indicating  short-crestedness  (Fig.  9a  and  c). 

Fig.  10  shows  the  quantitative  comparison  of  our  model  results  to  those  from  the  time- 
dependent  Boussinesq  model  of  Chen  et  al.  (1999)  along  the  central  cross-shore  transect. 
The  surface  elevation  from  two  model  simulations  agrees  reasonably  well  (Fig.  10a). 
Minor  differences  are  possibly  due  to  the  nonlinear  shallow  water  effects  (wave-wave 
interaction)  that  are  not  included  in  (8).  The  Boussinesq  model  surface  elevations  show 
sharp  crests  and  wide  troughs  in  the  surf  zone  (Fig.  10a),  indicating  the  importance  of  the 
nonlinearities  in  the  shallow  water  areas.  This  results  in  a  phase  offset  between  the  two 
modeled  water  surface  elevations  in  the  shallow  region.  The  accumulated  phase  offset 
increases  in  the  shoreward  direction.  Panchang  and  Demirbilek  (2001)  describe  the  finite 
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Fig.  9.  Model  results  for  the  rip  current  case  of  Chen  et  al.  (1999).  (a)  Surface  elevation;  (b)  amplitude;  (c)  phase. 


element  modeling  of  (8)  when  wave-wave  interactions  are  included;  however,  the 
combination  of  wave-wave  interaction  and  wave-current  interaction  has  not  yet  been 
reported.  Nevertheless,  we  have  found  that  incorporating  nonlinear  dispersion  relation  (5) 
can  partially  remedy  this  situation.  This  is  shown  in  Fig.  10  (bottom  panel),  where  the 
phase  offset  has  been  largely  eliminated.  Discrepancies  in  Fig.  10  may  also  be  attributed  to 
differences  in  the  implementation  of  wave  breaking  in  the  models.  The  amplitude 
distribution  predicted  by  the  present  model  (shown  as  a  dotted  curve  in  Fig.  10a)  identifies 
the  breaking  point  to  be  at  v  =  1 1 .2  m,  while  results  from  the  Boussinesq  model  show  wave 
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Fig.  10.  Simulations  of  the  rip  current  case  of  Chen  et  al.  (1999). 


heights  decreasing  only  after  x=\2m.  The  absence  of  benchmark  data  is  a  handicap 
preventing  us  from  further  verifying  the  breaking  criteria  used  the  two  models. 


6.  Wave  pattern  near  a  jettied  tidal  inlet 

The  effect  of  tidal  currents  on  wave  propagation  near  an  entrance  to  a  harbor  is  a 
problem  of  interest  in  coastal  engineering  practice  since  it  affects  navigation  safety,  inlet 
entrance  maintenance,  jetty  design,  harbor  performance,  etc.  Smith  et  al.  (1998)  conducted 
a  laboratory  study  to  investigate  the  wave-current  interaction  and  the  associated  wave 
breaking  in  the  presence  of  ebb  currents  in  an  idealized  entrance.  We  use  this  case  as  an 
example  to  develop  an  understanding  of  the  wave  patterns  that  may  be  expected  in  such 
configurations.  Since  the  hydraulic  model  measurements  for  waves,  currents,  and 
bathymetry  were  not  detailed  enough  for  a  systematic  model  validation,  a  comparison  can 
be  performed  only  in  a  general  sense.  However,  the  comparison  is  sufficiently  satisfactory 
and  the  numerical  results  are  sufficiently  informative  (in  fact,  they  shed  more  light  on  the 
phenomena  than  the  laboratory  data  do)  that  they  provide  the  motivation  to  supplement  the 
hydraulic  model  study  with  numerical  investigations  of  additional  conditions  of  interest  in 
this  harbor  engineering  problem. 

6.1.  Summary  of  hydraulic  model  data 

Fig.  1 1  (Smith  et  al.,  1998)  shows  the  overall  arrangement  of  the  physical  model  facility 
built  at  a  prototype  scale  of  1:50.  A  wave  maker  is  situated  24  m  offshore  at  the  head  of 
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Fig.  11.  Hydraulic  model  layout,  idealized  inlet  facility  (after  Smith  et  al.,  1998). 

the  model  basin  and  currents  can  be  generated  by  a  pumping  system.  Fig.  12  shows  the 
detailed  geometry  and  bottom  topography  in  the  vicinity  of  the  jettied  inlet,  where  model 
simulations  are  focused.  The  idealized  bathymetry  is  based  on  the  design  data  for  physical 
model  construction.  Most  of  the  domain  consists  of  a  beach  that  slopes  down  to  a  depth  of 
about  0.2  m,  followed  by  a  steep-slope  leading  to  the  deepwater  end,  where  water  depth  is 
constant.  Two  parallel  jetties  are  built  to  create  an  entrance  channel,  between  which  the 
depth  is  constant.  This  is  preceded  by  a  deep  trench  at  the  throat,  which  connects  the  backbay 
area  (harbor)  to  the  entrance  channel.  Experiments  were  limited  to  normally  incident  waves 
and  ebb  currents.  Wave  heights,  currents,  and  water  depths  were  measured  at  six  gauges 
(Fig.  12)  along  the  main  channel.  Laboratory  runs  were  made  using  narrow-band 
unidirectional  waves  with  desired  target  significant  wave  heights  and  peak  periods.  These 
are  simulated  numerically  using  a  monochromatic  representation  with  the  associated 
significant  wave  heights  and  peak  periods.  This  assumes  that  narrow-band  unidirectional 
seas  can  be  reasonably  approximated  by  the  equivalent  monochromatic  waves. 


Fig.  12.  Wave  model  domain,  detailed  bathymetry  (depth  in  m),  and  underlying  ebb  current.  Circles  indicate 
gauge  locations  numbered  from  one  at  offshore  to  six  near  the  inlet. 
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6.2.  Flow  field  simulation 

A  detailed  current  field  was  not  available  from  hydraulic  model  study  since 
measurements  were  made  only  at  the  locations  shown  in  Fig.  12.  A  synthetic  flow  field 
had  to  be  generated  numerically  to  approximate  the  actual  flow  conditions  in  the 
laboratory.  The  two-dimensional  circulation  model  DUCHESS,  developed  at  TU  Delft 
(Booij,  1989;  Jin  and  Kranenburg  1993)  was  applied  to  a  domain  extending  +  30  m  in  the 
longshore  direction  and  30  m  to  the  offshore  direction,  large  enough  to  minimize  any 
potential  adverse  lateral  open  boundary  effects.  The  backbay  area  is  not  of  interest  for  this 
study  and  was  eliminated  by  placing  a  water  boundary  at  the  inlet  throat  (y  =  —  4  m).  The 
volume  flux  along  this  boundary  with  certain  distribution  was  prescribed  to  drive  the 
model  for  desired  ebb-flow  field.  All  offshore  model  boundaries  are  treated  as  flow¬ 
through  (non-reflecting)  using  a  radiation  boundary  condition.  To  prevent  instabilities  in 
the  model  solution,  the  volume  flux  was  ramped  up  gradually  from  zero  to  a  specified 
value,  and  was  maintained  constant,  thereafter,  until  the  current  field  in  the  entire 
modeling  domain  reached  equilibrium.  Input  volume  flux  (which  is  unavailable  from  the 
hydraulic  model  data),  bottom  friction,  and  eddy  viscosity  were  used  as  calibration 
coefficients  in  order  to  obtain  the  best  fit  between  simulated  currents  and  the  data  observed 
at  the  gauges. 

The  resulting  jet- like  ebb  current  field  from  model  simulation  is  shown  in  Fig.  12.  (For 
clarity,  only  a  part  of  the  flow  model  domain  is  shown).  Detailed  velocity  profiles  are 
provided  in  Fig.  13a  along  four  transects  across  the  mainstream  of  the  ebb  flow.  The 
velocity  distributions  show  a  common  structure  of  double-peak  shear  currents,  with  the 
shear  structure  intensifying  from  the  offshore  to  a  maximum  near  the  middle  of  the  throat, 
and  then  decreasing  as  a  result  of  the  specified  flowfield  being  constant.  Fig.  13b  shows  the 
variation  of  the  current  velocity  along  the  centerline.  A  rapid  change  in  velocity  at  near 
y  =  0  is  due  to  the  flow  transition  over  a  steep-slope  from  deep  trench  in  the  inlet  throat  to  a 
shallow  flat  area  between  the  jetties.  Also  in  Fig.  13b,  the  simulated  currents  are  compared 
to  the  measurements  of  two  laboratory  runs  (Runs  3  and  4,  described  later)  with  the  same 
target  flow  conditions,  showing  apparently  good  agreement  in  both  cases.  However,  it  is 
still  uncertain  how  well  the  simulated  current  field  represents  the  actual  flow  in  the 
laboratory  study.  In  particular,  the  limited  data  prevent  us  from  confirming  the  lateral 
(shearing)  structure  of  the  ebb  flow  seen  in  Fig.  13a. 


Fig.  13.  Simulated  current  velocity  profiles,  (a)  Longshore  transects;  (b)  cross-shore  transect,  compared  to  lab 
data. 
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6.3.  Simulation  of  normal  wave  incidence 


We  investigated  two  cases,  without  currents  and  with  ebb  currents,  each  of  which  consists 
of  a  large-wave  case  and  a  small-wave  case.  The  target  significant  incident  wave  amplitude  is 
2.75  cm  for  large-waves  and  1.83  cm  for  small-waves  in  the  laboratory  study.  The  wave 
period  is  1 .4  s  for  all  cases.  The  computational  domain  for  wave  simulations,  shown  in  Fig.  12, 
consists  of  165,371  triangular  elements  and  83,581  nodes  (representing  a  resolution  of  20 
nodes  per  wavelength),  and  the  radius  of  the  semicircular  boundary  is  13.6  m.  All  shoreline 
and  jetty  boundaries  are  considered  non-reflecting  and  assigned  the  new  boundary  condition 
(15)  for  absorbing  outgoing  waves.  The  advantage  of  the  direction-dependent  boundary 
condition  described  in  Section  2  can  be  demonstrated  in  this  application,  since  the  incident 
wave  angles  along  much  of  the  jetty  boundaries  are  likely  to  be  large.  This  is  demonstrated  in 
Fig.  14  for  the  case  of  normal  wave  incidence  and  no  ambient  current.  Fig.  14b  shows  that 
spurious  reflections  seen  in  Fig.  14a  beyond  the  tips  of  the  jetties  have  disappeared,  and  waves 
can  follow  the  sidewalls  of  the  jetties  smoothly  as  they  propagate  into  the  inlet.  Fig.  15 
compares  these  results  with  hydraulic  model  measurements  and  indicates  that  model 
estimates  with  the  new  boundary  condition  have  a  somewhat  better  match  to  the  data. 
Obviously,  with  the  boundary  condition  (15)  that  accommodates  outgoing  wave  directionality 
in  a  nonlinear  manner,  the  model  estimates  are  much  improved. 

Model  results  for  the  cases  of  no  currents  are  presented  and  compared  to  the  data  in 
Fig.  16.  Runs  1  and  2  correspond,  respectively,  to  large-  and  small-wave  incidence. 
Because  of  the  uncertainties  in  the  hydraulic  model  associated  with  establishing  the 
targeted  (as  opposed  to  the  actual)  incident  wave  heights,  different  incident  values  were 
tried  until  the  best  match  was  achieved  at  the  offshore-most  gauge  (Gauge  1).  Fig.  16 
shows  that  model  predictions  overestimate  the  data  when  the  target  significant  wave 
heights  are  used  as  the  incident  wave  conditions.  With  smaller  amplitudes  At  =  2.5  and 
1.65  cm  for  Runs  1  and  2,  the  best  results  are  obtained  relative  to  the  data  at  all  gauges 
except  for  Gauge  2.  Note  that  data  at  Gauge  2  are  consistently  lower  than  model  estimates 
for  all  conditions  (Fig.  18,  see  Section  6.4  for  further  discussion). 


Fig.  14.  Simulated  wave  amplitudes,  jettied  inlet  without  current,  (a)  Conventional  boundary  condition;  (b) 
improved  boundary  condition. 
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Fig.  15.  Model  results  of  wave  amplitudes  at  central  transect,  compared  to  lab  data. 

For  ebb  currents,  model  estimates  were  made  using  adjusted  wave  amplitudes  of  2.5  cm 
(Run  3)  and  1.65  cm  (Run  4).  Based  on  the  results  for  Run  3  (Fig.  17),  the  model  captures 
the  bending  effect  inflicted  on  the  crests  by  the  opposing  currents  (Fig.  17a);  the  phase 
pattern  in  the  middle  of  the  entrance  channel  indicates  wave  crossing  and  focusing  (as  in 
the  case  of  shoals,  Berkhoff  et  al.,  1982).  Waves  heights  are  amplified  not  only  by  the 
opposing  currents,  but  also  by  the  increasing  focusing  due  to  current  shear,  which  creates  a 
high-wave  zone  in  the  middle  of  the  domain  (Fig.  17b).  The  width  of  the  high- wave  zone 
decreases  along  the  centerline  towards  the  middle  of  the  entrance  channel,  where  waves 
cross,  leaving  large  shadow  areas  near  the  jetty  walls.  Fig.  18  shows  that,  for  both  large- 
and  small- wave  incidence  (Runs  3  and  4),  following  the  centerline,  waves  increase  rapidly 
in  amplitude  and  become  large  enough  to  break.  The  first  wave  breaking  occurs  further 
offshore  (y  =  6.5  m)  for  the  large- wave  case,  compared  to  y  =  4.5  m  for  both  the  small- 
wave  case  and  the  no  current  case  (Fig.  16a).  Continued  strong  focusing  by  the  underlying 
currents  in  the  deep,  flat,  and  narrow  channel  between  the  jetties  results  in  renewed  growth 
until  a  second  wave  breaking  occurs  at  y  =  2.4  m. 

Model  results  for  Runs  3  and  4  provide  fair  agreement  with  the  laboratory  data 
(Fig.  18),  although  they  are  not  as  good  as  that  without  currents;  the  predicted  wave 
amplitudes  are  generally  higher  than  the  data.  As  noted  above,  the  uncertainty  of  the  actual 
current  field  used  in  the  model  compared  to  the  laboratory  condition  is  perhaps  the  main 
reason  for  the  discrepancies.  Also,  we  were  unable  to  explain  why  Gauge  2  measurements 
for  all  test  conditions  were  consistently  lower  than  model  estimates.  If  instrumentation 


Fig.  16.  Wave  amplitudes  along  the  centerline,  jettied  inlet  case,  normally  incident  waves  without  current. 
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Fig.  17.  Model  results  for  a  jettied  inlet  case,  with  ebb  current,  (a)  Phase;  (b)  amplitude. 


errors  were  involved  at  this  gauge,  we  would  have  lost  important  information  about  the 
onset  of  wave  breaking.  Yet,  as  Fig.  17  reveals,  an  intricate  wave  pattern  may  be  expected 
near  harbor  entrances  with  jetties;  it  also  suggests,  however,  that  many  more 
measurements  would  be  needed  to  identify  any  meaningful  features  of  the  wave-current 
interaction  in  the  hydraulic  model. 


6.4.  Oblique  wave  incidence 

Since  there  are  no  data  available  for  oblique  wave  incidence  from  the  hydraulic  model 
study,  we  investigate  the  effect  of  ebb  and  flood  currents  on  wave  transformation.  Fig.  19 
shows  numerical  results  for  two  simulations,  i.e.  for  45°  incidence  (A,  =  2.5  cm) 
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Fig.  18.  Wave  amplitudes  along  the  centerline,  jettied  inlet  case;  normal  waves  with  ebb  currents. 
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Fig.  19.  Simulated  amplitudes  for  oblique  wave  incidences,  jettied  inlet  case,  (a)  With  ebb  current;  (b)  with  flood 
current. 

encountering  an  ebb  current  and  a  flood  current.  The  representing  flood  current  is  assumed 
to  be  simply  the  reverse  of  the  ebb  current  for  the  purposes  of  this  case  study. 

Oblique  waves  are  refracted  to  varying  degrees  when  they  come  across  the  ebb  flow 
that,  as  shown  in  Figs.  12  and  13,  has  a  structure  with  changing  shear  intensity.  As  shown 
in  Fig.  19a,  the  double-peak  shear  current  causes  wave  focusing  and  diverging,  which 
generates  two  elongated  low-wave  and  two  high- wave  zones  in  the  central  area.  Fig.  20b 
clearly  indicates  this  large  fluctuation  in  the  wave  amplitudes  across  the  ebb  stream,  with 
two  peaks  corresponding  to  wave  focusing  due  to  increasing  current  shear.  Compared  to 
the  case  of  normal  wave  incidence  (Fig.  20a),  peak  waves  in  the  flow  zone  are  much 
weaker  (also  see  Fig.  21). 

With  flood  currents,  the  strong  current  shear  near  the  left  edge  of  the  flood  stream  does 
not  allow  oblique  waves  to  fully  penetrate  into  the  flow  zone  (McKee,  1974). 
Consequently,  part  of  the  wave  energy  is  reflected  toward  the  left  (Fig.  19b),  causing 
large  fluctuations  of  wave  amplitude  on  the  left  side  of  the  mainstream  (Fig.  20b).  The  fact 
that  similar  effects  were  obtained  in  the  vortex  ring  test  lends  credence  to  these  results.  The 
transmitted  part  of  waves  is  subject  to  refraction  by  the  flood  current,  and,  as  in  the  case 
with  ebb  currents,  results  in  large  variation  in  wave  amplitudes  inside  the  flow  zone. 
For  flood  currents,  the  wave  amplitudes  are  much  lower  in  the  flow  zone  outside  the  tip  of 
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Fig.  20.  Simulated  wave  amplitudes  along  two  lateral  transects,  jettied  inlet  case  (At  =  0.025  m). 
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Fig.  21.  Simulated  wave  amplitudes  along  the  centerline,  jettied  inlet  case  (A,-  =  0.025  m). 


jetties  because  of  wave  reflection  at  the  left  edge  of  the  mainstream  (Figs.  19b  and  21)  and 
also  because  the  waves  are  following  the  currents.  In  spite  of  this,  more  wave  energy  is 
able  to  penetrate  into  the  jetty  channel  and  the  inlet  compared  to  the  ebb  current  case, 
either  with  normal  or  oblique  incidence  (Fig.  21).  This  is  because  flood  currents  have 
turned  the  waves  more  toward  the  inlet,  and  that  lower  wave  amplitudes  have  resulted  in 
less  wave  energy  dissipation  (breaking)  in  the  channel  and  inlet  area  between  jetties. 


7.  Concluding  remarks 

We  have  provided  a  comprehensive  treatment  of  the  various  issues  associated  with 
finite-element  modeling  of  the  elliptic  wave-current  interaction  equation.  Methods  to 
effectively  handle  various  nonlinear  aspects  of  the  solution  are  developed.  Implementation 
of  an  improved  boundary  absorption  condition  that  accounts  for  the  angle  of  wave 
approach  and  the  local  bathymetric  variation  has  largely  eliminated  the  spurious 
oscillations  commonly  seen  in  the  solutions  of  similar  models.  The  model  developed  on 
the  basis  of  the  methods  detailed  here  was  applied  to  three  wave-current  interaction  cases. 
The  results  were  found  to  be  consistent  with  those  of  other  published  models. 

Relative  to  the  frequently  used  parabolic  approximation  models  (e.g.  Yoon  and  Liu, 
1989),  the  present  results  are  nearly  identical  in  some  cases.  In  other  cases,  the  limitations 
of  the  parabolic  approximations  are  illustrated,  such  as  the  effect  of  dropping  the  second- 
order  derivative  when  the  wave  height  changes  rapidly  (as  in  the  rip-current  test-case). 
Also,  the  reflective  pattern  seen  in  the  elliptic  model  results  as  a  consequence  of  a  highly 
oblique  angle  at  which  waves  approach  the  current  field  permits  one  to  speculate  that 
parabolic  approximations  may  be  inadequate  in  the  case  of  strong  velocity  shears,  where 
large-angle  wave  refraction,  diffraction,  and  even  reflection  (due  to  wave  caustics)  may  be 
expected.  Relative  to  Boussinesq  models,  which  have  traditionally  been  applicable  only  to 
shallow  water  waves1  the  present  approach  suffers,  in  principle,  from  the  disadvantage  of 
inadequate  representation  of  nonlinear  shallow- water  effects.  However,  the  results 
obtained  here  show  only  minor  discrepancies,  and  even  those  are  somewhat  diminished 


The  applicability  of  Boussinesq  models  is  now  being  extended  to  intermediate  and  deep  water  depths. 
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through  the  use  of  the  nonlinear  dispersion  relation  (5).  Also,  a  strategy  for  including 
wave-wave  interactions  in  elliptic  finite  element  models  and  some  preliminary  results 
have  been  described  in  Panchang  and  Demirbilek  (2001).  Since  multiple  nonlinear 
mechanisms  have  been  successfully  handled  in  the  work  reported  here,  incorporating  the 
additional  mechanism  of  wave-wave  interaction  may  not  be  problematic.  Further,  a  finite- 
element  formulation  that  can  better  represent  complex  geometries  (such  as  the  case  of 
jettied  inlets)  is  more  conducive  to  the  time-independent  elliptic  formulation  (8)  than  to 
the  time-dependent  and  highly  nonlinear  Boussinesq  models. 

The  satisfactory  solutions  obtained  for  the  tests  described  in  Section  5  prompted  the 
investigation  of  wave-current  interactions  near  a  jettied  tidal  inlet,  a  problem  of 
considerable  engineering  importance  (e.g.  Lin  and  Demirbilek,  in  press;  Smith  et  al., 
1998).  The  wave  patterns  are  seen  to  be  quite  intricate.  Ebb  currents  lead  to  considerable 
breaking  outside  the  entrance  channel,  but  renewed  wave  growth  occurs  due  to  focusing  by 
the  currents  before  the  waves  break  again.  For  flood  tides,  the  wave  heights  are  smaller 
outside  the  jetty  channel,  as  a  consequence  of  the  non-opposing  nature  of  the  waves  and 
currents,  but  they  are  larger  than  in  the  ebb  current  case  inside  the  channel  because  waves 
do  not  break.  Although  the  case  investigated  here  dealt  with  the  layout  constructed  by 
Smith  et  al.  (1998),  the  overall  results  presented  here  have  set  the  stage  for  tackling 
practical  engineering  problems. 

The  solution  methods  proposed  here  require  only  moderate  computer  power.  The 
examples  illustrated  show  that  a  typical  application  involving  monochromatic  wave 
transformation  would  take  about  15  min  of  run  time  on  a  modem  PC. 
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